{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "chicken-gates",
   "metadata": {},
   "source": [
    "# Example 14.5 Pressure Relief"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "illegal-flight",
   "metadata": {},
   "source": [
    "A 2000 liter vessel contains 1000 kg of a 10 wt% water and 90 wt% propanediol-1,2 mixture, stored at its own vapor pressure, and protected by a rupture disk which is set to activate at 3 bar. A heat flux of 1 MW is applied to the fluid, which will cause a relief event. Determine the mass flow rate which will exit the rupture disk.\n",
    "\n",
    "\n",
    "Use the NRTL model with A12 = 0.1078, A21=-.2811, B12=62.0818, B21=-1.4101, alpha12=0.3, and an ideal gas model for the vapor phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "imported-shield",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17378.254767083745"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from thermo import ChemicalConstantsPackage, NRTL, IdealGas, GibbsExcessLiquid, FlashVL\n",
    "from chemicals import ws_to_zs, mixing_simple, rho_to_Vm\n",
    "from scipy.constants import bar\n",
    "# Propylene glycol is another more common synonym than propanediol-1,2\n",
    "constants, correlations = ChemicalConstantsPackage.from_IDs(['water', 'Propylene glycol'])\n",
    "from scipy.constants import calorie, R\n",
    "zs = ws_to_zs(ws=[.1, .9], MWs=constants.MWs)\n",
    "\n",
    "\n",
    "P = 1*bar\n",
    "T = 300 # Initial guess\n",
    "V_vessel = 2.0 # m^2\n",
    "m = 1000 # kg in the tank\n",
    "\n",
    "# Calculate the average molecular weight\n",
    "MW_avg = mixing_simple(zs, constants.MWs) # g/mol\n",
    "# Compute the overall density for the contents of the tank\n",
    "bulk_density = m/V_vessel # kg/m^3\n",
    "# Compute the molar volume of the bulk solution\n",
    "Vm_avg = rho_to_Vm(rho=bulk_density, MW=MW_avg) # m^3/mol\n",
    "moles = V_vessel/Vm_avg # Compute the total number of moles of substance in the tank\n",
    "moles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "coastal-present",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure the liquid phase\n",
    "# As this system is low-pressure, the Poynting correction factor is ignored\n",
    "# Since the gas phase is set to the idea gas, it makes no sense to use\n",
    "# Phi correction factors for gas-nonideality in the liquid phase either\n",
    "tau_as = [[0, 0.1078/R*calorie], [-.2811/R*calorie, 0]]\n",
    "tau_bs = [[0, 62.0818/R*calorie], [-1.4101/R*calorie, 0]]\n",
    "alphaC =  [[0, 0.3],[.3, 0]]\n",
    "GE = NRTL(T=T, xs=zs, tau_as=tau_as, tau_bs=tau_bs)\n",
    "\n",
    "liquid = GibbsExcessLiquid(VaporPressures=correlations.VaporPressures,\n",
    "                           HeatCapacityGases=correlations.HeatCapacityGases,\n",
    "                           VolumeLiquids=correlations.VolumeLiquids,\n",
    "                           GibbsExcessModel=GE,\n",
    "                           equilibrium_basis='Psat', caloric_basis='Psat',\n",
    "                           T=T, P=P, zs=zs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "planned-wyoming",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure the ideal gas\n",
    "gas = IdealGas(T=T, P=P, zs=zs, HeatCapacityGases=correlations.HeatCapacityGases)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "soviet-farming",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "405.05615108404197"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import brenth\n",
    "\n",
    "flasher = FlashVL(constants, correlations, liquid=liquid, gas=gas)\n",
    "# To find the temperature at which the rupture disk will break open, we have two specs\n",
    "# 3 bar and the molar volume of the bulk solution calculated earlier, `Vm_avg`.\n",
    "# Iteration on PT flashes is one solution.\n",
    "\n",
    "def V_error(T):\n",
    "    PT = flasher.flash(T=T, P=P, zs=zs)\n",
    "    return PT.V() - Vm_avg\n",
    "\n",
    "T_release = brenth(V_error, 300, 600, xtol=1e-6)\n",
    "T_release"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "passing-arizona",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "131.906151084042"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T_release-273.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "descending-queue",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['HEOS_FIT', 'DIPPR_PERRY_8E']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[i.method for i in correlations.VolumeLiquids]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "resistant-chart",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}